Excited states dynamics in time-dependent density functional theory: high-field 
molecular dissociation and harmonic generation. 
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We present a theoretical description of femtosecond laser induced dynamics of the hydrogen 
molecule and of singly ionised sodium dimers, based on a real-space, real-time, implementation of 
time- dependent density functional theory (TDDFT) . High harmonic generation, Coulomb explosion 
and laser induced photo-dissociation are observed. The scheme also describes non-adiabatic effects, 
such as the appearance of even harmonics for homopolar but isotopically asymmetric dimers, even 
if the ions were treated classically. This TDDFT-based method is reliable, scalable, and extensible 
to other phenomena such as photoisomerization, molecular transport and chemical reactivity. 
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It is now possible to study electron and molecular dy- 
namics in real time using various experimental techniques 
employing intense ultra-short laser sources |Q . Some ex- 
amples of such investigations include X-ray photoelec- 
tron spectroscopy of molecules pump-probe ionisa- 
tion measurements || , production of high harmonics as 
a source of soft X-rays Q, the measurement of electron- 
phonon interactions in thin films B, and the estimation 
of the onset of Coulomb screening || . A technologically 
important and very active field of research is the appli- 
cation of ultra-short laser pulses to induce, control and 
monitor chemical reactions [j?], ^, Whenever the in- 
tensity of the laser field is comparable to the molecular 
electronic fields, perturbative expansions break down and 
new processes appear, which are not fully understood 
from a microscopical point of view pp| . A practical and 
accurate computational framework to descibe excited- 
state electron-ion dynamics is therefore still needed. 

Not surprisingly, the smallest systems have attracted 
particular attention from both experimentalists and theo- 
reticians, as a bench-horse to improve our understanding 
of electron dynamics at the femtosecond scale |ll], [l2| . 
However, the methods used in these calculations can not 
be easily extended to larger and more realistic systems. 
The exact quantum mechanical solution of a 3D system 
of more than three particles is certainly not feasible with 
state-of-the-art computers. ID models are much easier to 
handle, but they can not really be used as predictive tools 
for problems involving the interaction of lasers with large 
clusters or solid-state systems of technological relevance. 

To tackle such a problem, time-dependent density 
functional theory (TDDFT) appears as a valuable 
tool. Even with the simplest approximation to the 
exchange-correlation potential, the adiabatic local den- 
sity approximation (ALDA), one obtains a very good 
compromise between computational ease and accuracy 
Uj . TDDFT can certainly be applied to large systems 
in non perturbative regimes, while providing a consistent 
treatment of electron correlation. It has been well tested 



in the study of electron excitations, like the optical ab- 
sorption spectra in the linear regime [2l[] . Although 
almost all applications of TDDFT in the field of laser 
physics have only involved electronic dynamics, recent 
attempts have also been made at describing the coupled 
nuclear and electronic motion in laser fields |E| , account- 
ing for the nuclear motion classically. A full quantum 
mechanical treatment of the system could in principle be 
done within a multi-component TDDFT, although it has 
not been tried for more than three particles p3fl . How- 
ever, since many vibrational quanta are coherently ex- 
cited, there is good motivation for the classical treatment 
of the nuclei. The purpose of this work is to illustrate a 
general method to study many-electron systems subject 
to strong laser fields. It is based on the quantum mechan- 
ical propagation of the electronic wave packet - described 
within TDDFT - combined with classical motion of the 
nuclei. As the laser field populates the excited Born- 
Oppenheimer surfaces, this scheme includes diabatic ef- 
fects, while maintaining a good scaling with the size of 
the system. As an illustration we focused on one and two 
electron dimers, namely NaJ and the hydrogen molecule. 

The equations of motion may be derived from the La- 
grangian: 
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where Edft is the usual Kohn-Sham density functional, 
depending on the electron orbitals {<fi} and the nuclear 
coordinates {R}, and £(t) is the time-dependent electric 
field from the laser pulse. Variation of the Lagrangian 
then yields Newton's equations for the nuclear coordi- 
nates, 
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and the usual TDDFT equations for the orbital variables 
|L8||. We solved these equations in real time, following the 
method of Yabana and Bertsch BO], using a real-space 
grid representation of the orbitals p4 pH| . This scheme 
has the advantage that the Kohn-Sham Hamiltonian is 
a very sparse matrix. The forces in Eq. (||) are calcu- 
lated with the help of a generalised Hellmann-Feynman 
theorem, 
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where Hks = SEdft / 8n(r, t) is the Kohn-Sham Hamil- 
tonian. For numerical reasons we represent the electron- 
ion interaction by norm-conserving non-local Troullier- 
Martins pseudopotentials |2(|. 

We have studied two different classes of time- 
dependent problems: photofragmentation and high har- 
monic generation. We now discuss the first case, the Naj" 
dimer in a femtosecond laser field. It is a good test, since 
it has been exhaustively studied using a diverseness of 
approaches (|, [l4|, |27|. In particular, a recent experi- 
ment J8| focused on the photofragmentation of NaJ in 
intense femtosecond laser fields. Using a pump-probe 
technique, the authors discovered that NaJ dissociated 
in four different channels, ranging from simple field ioni- 
sation followed by Coulomb explosion, to photodissocia- 
tion on light-induced potentials. For these calculations, 
we used an uniform grid spacing of 0.3 A, and the sys- 
tem was confined to a sphere of radius 10 A. In these 
one-electron calculations, we omitted the core-valence 
exchange-correlation. In this case the total electronic 
energy is given by the Kohn-Sham eigenvalues, and they 
can be used to compute the adiabatic potential energy 
surfaces shown in Fig. HI. The two lowest single-photon 
transitions from the 1^^" ground state are at 2.5 eV 
(to the 1 2 £+ state) and at 3.2 eV (to the 1 2 IT U state). 
The latter is achieved by a laser polarized perpendicu- 
larly to the internuclear axis. These energies accord well 
with the observed single-photon transitions [^8) . For the 
time-dependent calculations, we start with the dimer in 
its ground state, which is then propagated with a modi- 
fied Krank-Nicholson scheme . The time step for the 
time integration was 0.005 h eV^ 1 ks 0.003 fs. A simple 
check on the implementation of the time evolution op- 
erator consists of calculating the linear photo-absorption 
spectrum, using a weak ^-function external field, as in 
Ref. |^0| . Almost all the spectral weight is concentrated 
in two peaks (see inset in Fig. |l|), which are at energies 
corresponding exactly to the vertical transitions between 
energy surfaces. We note that this exact correspondance 
is only obtained for one-electron systems: in general the 
TDDFT spectra will have shifts from the energy surfaces 
determined by the Kohn-Sham eigenvalues. 

Next we examine the evolution of the dimer under 
high-field excitation. We consider external fields of the 
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FIG. 1: Adiabatic energy surfaces of the NaJ dimer, as 
obtained by our three-dimensional real-space code. Similar 
results can be found in Ref. [p9|| . In the inset, the computed 
photoabsorption cross-section of the same molecule. 
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where Iq is the maximum intensity of the pulse, e is the 
polarization vector, and r is the pulse duration, taken as 
r = 80 fs. As a first case, we examine the effect of excita- 
tion at the lower resonant frequency, lo = 2.5 eV. In Fig. 
2(a), we present a series of runs at different intensities, 
ranging from weak (10 10 W/cm 2 ) to moderate (2.1 x 10 12 
W/cm 2 ). Since the l 2 Sj surface is anti-bonding, exci- 
tation at this resonant frequency should lead to dissoci- 
ation, even at moderate intensities. This is indeed con- 
firmed by our calculations. The upper panel depicts the 
internuclear separation of the dimer, which exhibits an 
acceleration during the laser pulse and a nearly constant 
velocity expansion thereafter. Clearly, the dimer disso- 
ciates at all field levels that we applied. To examine 
the ionisation of the dimer, we assumed that any den- 
sity reaching the edges of the simulation box corresponds 
to unbound electrons. By absorbing this density at the 
boundaries, we can thus define the ionisation probability 
as I(t) = 1 — N(t), where N(t) is the charge that remains 
inside the simulation box at time t. The lower panel of 
Fig. 2(a) shows A as a function of t. We see that there is 
practically no ionisation for the lower fields, and only a 
20% ionization probability for the 2.1 x 10 12 W/cm 2 field. 
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FIG. 2: Evolution of internuclear distance (top panel) and electronic charge in simulation region (bottom panel) for the NaJ 
molecule. The dimers are excited with laser pulses of 2.5, 3.2 and 1.57 eV in columns (a), (b) and (c) respectively. 



Thus, in this range of intensities, the laser dissociates the 
dimer without ionising it. 

We next consider the excitation at the upper resonance 
frequency, u) = 3.2 eV, corresponding to an electric field 
perpendicular to the dimer axis. Since the l 2 n„ surface 
is bonding, [see Fig. || (b)] no dissociation is expected 
unless the Coulomb explosion channel is opened through 
ionisation. We see that the dimer remains bound over 
the entire range of intensities that produced dissociation 
at the lower resonant frequency. 

Finally, we also performed simulations at the non- 
resonant frequency lu = 1.57 eV, the one used in Ref. 
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Harmonic Order 

FIG. 3: Harmonic spectra of HD (left panels) and H2 (right 
panels) . The nuclear masses used in the calculation are m' = 
am, being m the real mass. In this way, top plots were made 
using for the nuclear masses their real values whereas bottom 
plots were made using a hundredth of their real values. 



H . Fig. H (c) shows how dissociation now occurs only at 
much higher intensities, and it is mainly due to ionisation 
(almost absent in the resonant calculations for the range 
of intensities used): since Na 2 , 4 " has no bonding states, 
ionisation is followed by Coulomb explosion. 

Another process in which the nuclear motion may play 
an important role is high harmonic generation. Even har- 
monics may be created by irradiating HD with an intense 
laser pulse, but not by irradiating H2: even harmonic 
generation is forbidden for a centrosymmetric molecule. 
In an adiabatic treatment of the nuclear coordinates, the 
nuclear masses play no role and the even harmonics can 
not appear. This is no longer the case if non-adiabatic 
effects are taken into account, for the different masses of 
H and D break the symmetry. Kreibich et al. |l7| studied 
this process in a 1-D model with a full quantum mechan- 
ical treatment of the nuclear motion, finding strong even 
harmonics at high harmonic number. To discern whether 
the classical treatment of nuclear motion also produces 
these harmonics, we studied the same ID problem within 
our framework. As in Ref. |l7| , we took the laser field 
to have a frequency of 1.6 eV, and an intensity that rises 
linearly to 10 14 W/cm 2 over an interval of 10 optical cy- 
cles, and is held constant thereafter. We then calculated 
the spectral intensity of the generated harmonics, H(lu): 
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We find that the classical treatment does indeed produce 
even harmonics, but much smaller than the quantum 
treatment. The results are shown in Fig. 0. The top 
left panel depicts the harmonic spectrum for HD, and 
only odd harmonics are apparent. However, it may be 
proved that the HD Hamiltonian already violates cen- 
trosymmetry within our classical treatment, through a 
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term of the form: 

where P(t) = | (P# (t) — Po(t)) is the relative time- 
dependent nuclear momentum and p$ are the electronic 
momentum operators . Its effect can be enhanced by 
decreasing the nuclear masses. In the bottom left panel, 
the H and D masses have been decreased by a factor 100, 
and then the second- and fourth-order harmonics become 
visible. As a qualitative check of the numerics, we also 
show the same graphs for H2 , in which no even harmonics 
can occur. 

Thus we see that on qualitative level the non-adiabatic 
dynamics generating even harmonics are obtained with 
the classical treatment of the nuclear coordinates. How- 
ever, the quantum treatment may be needed for a quanti- 
tative result. By describing the nuclei quantum mechan- 
ically, the ground state violates centrosymmetry and the 
even harmonics can be generated even if the nuclear mo- 
tion is frozen. In contrast, in the classical treatment the 
ground state is symmetric and the symmetry violation 
only builds up as the nuclei move. 

In summary, we have examined the computational fea- 
sibility of including nuclear dynamics in time-dependent 



density functional theory using a pseudopotential code 
to study the femtosecond laser induced dynamics of 
sodium dimers. Using this approach for treating the NaJ 
dimer, we were able to distinguish different phodissocia- 
tion regimes, ranging from dissociation on light induced 
potentials, to field ionisation followed by Coulomb explo- 
sion. Electronic and ionic degrees of freedom are thus 
coupled, so that one can observe the electron-phonon 
transfer of energy. We also found, with another exam- 
ple, that non-adiabatic effects are present in the general 
treatment based on Eq. (1). One of the major attrac- 
tivenesses of this method resides in its reasonable scaling 
behaviour when applied to larger systems. We thus ex- 
pect to be able to tackle problems like photoisomerization 
or even photochemical reactivity in systems of dozens of 
atoms in the near future. 
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